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Abstract 

We study how the oscillations of the neutrinos affect their thermalization process during the 
reheating period with temperature O(l) MeV in the early universe. We follow the evolution of the 
neutrino density matrices and investigate how the predictions of big bang nucleosynthesis vary with 
the reheating temperature. For the reheating temperature of several MeV, we find that including 
the oscillations makes different predictions, especially for 4 He abundance. Also, the effects on the 
lower bound of the reheating temperature from cosmological observations are discussed. 
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I. INTRODUCTION 



The standard big bang model assumes that the universe was once dominated by thermal 
radiation composed of photons, electrons, neutrinos, and their anti-particles. It is one of the 
main issues in theories beyond the standard cosmology where these particles came from, or 
equivalently, what reheated the universe. The reheating temperature, at which the universe 
becomes radiation-dominated, is therefore a very important parameter that discriminates 
among many scenarios on the thermal history of the universe. In the following we define the 
reheating temperature as that of the latest reheating process, if the universe experienced 
several reheating stages. 

Recent observations of the cosmic microwave background radiation (CMB) Q has 
strongly suggested that the universe underwent inflation at an early stage. After infla- 
tion ended, the universe was dominated by the oscillation energy of the inflaton until it 
decayed and reheated the universe. The upper limit on the reheating temperature was ob- 
tained P| by constraining the relic abundance of the gravitinos, the superpartner of the 
graviton, which are inevitably present in the supersymmetric (SUSY) framework. Here we 
are interested in the relatively low reheating temperature, especially in the MeV range, and 
would like to put a lower limit on the reheating temperature. 

The MeV-scale reheating is actually ubiquitous in theories beyond the standard cosmol- 
ogy. In the framework of the SUSY and superstring theories, there are many particles with 
very long lifetime, e.g., the moduli and the gravitinos mentioned above, since their interac- 
tion is so weak, typically suppressed by the Planck scale. These long-lived massive particles 
might have dominated over the radiation from the inflaton decay. If the masses of these par- 
ticles are heavy enough, they decay and reheat the universe again just before the big bang 
nucleosynthesis (BBN) starts. Otherwise they often cause cosmological disaster known as 
"cosmological moduli problem" 3,0,0] and "gravitino problem" jJO, tIIs!- The simplest 
solution of these problems is to dilute the unwanted relics by producing large entropy at a 
later time [3, Q|. In either case, the reheating temperature is very low and typically around 
MeV. 

Another example that prefers the low reheating temperature is the curvaton scenario [h| 
in which the curvaton field dominates the universe and its isocurvature fluctuation is trans- 
formed into an adiabatic one. Furthermore, in the Affleck-Dine mechanism |l2j responsible 



2 



for the orig in of the baryon asymmetry, it is known that non-topological solitons such as 
Q-balls jl3[ are generally created. Since the decay process of the Q-balls is geometrically 



suppressed, they might dominate t 
studied in many different scenarios 



re universe, and such possibility has been extensively 



What if the reheating temperature is several MeV? In contrast to electrons that are always 
(at least until the temperature drops below a few eV) in thermal contact with photons via 
electromagnetic forces, neutrinos interact with electrons and themselves only through the 
weak interaction. The decoupling temperature of the neutrinos should be around 3 MeV for 
the electron neutrinos and 5 MeV for the muon and tau neutrinos, respectively 3, [3, Q • 
The difference comes from the fact that the electron neutrinos have additional charged 
current interaction with electrons. Therefore the neutrinos might not be fully thermalized 
if the reheating temperature is in the MeV range. If this is the case, the expansion rate 



of the universe becomes sma 
angular power spectrum [21 



ler, which affects the light element abundances and the CMB 
In particular it has been widely believed or taken 



for granted that the predicted abundance of 4 He decreases as the reheating temperature 
drops below a few MeV. This is because the smaller expansion rate delays the decoupling of 
the neutron-proton transformation, decreasing the neutron-to-proton ratio at the beginning 
of the BBN. Since almost all the neutrons are absorbed in the 4 He nuclei, such a naive 
reasoning can explain the dependence of the 4 He abundance on the reheating temperature. 
In this paper, however, we will see that this widespread picture is drastically changed if we 
take account of the neutrino oscillations. 

Recent neutrino oscillation experiments j3, H] have determined the mass differences and 
mixing angles with high precision and established that mixing angles are large. The crucial 
point is that a flavor eigenstate transforms itself into another one. Therefore we must take 
special care to calculate neutrino distribution functions and the resultant effective number 

22 , 3- S- 13 ) it i s useful t° follow the evolution of the 



of neutrinos. As pointed out in Refs. 



neutrino density matrices when flavor mixings are present. Here we will solve momentum 
dependent Boltzmann equations for the neutrino density matrices. We will see that the 
predicted abundance of 4 He is drastically changed, while the effective number of neutrinos 
does not change much. To put it simply, the reason for this is that the number density of the 
electron neutrinos is decreased due to flavor mixings, which makes the freeze-out temperature 
of the neutrons higher; this effect cancels and even overcomes that of the decrease in the 
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expansion rate. Thus MeV-scale reheating scenario is one of the examples in which the 
neutrino oscillations play a essential role. 

The outline of this paper is the following. In the next section we formulate neutrino 
thermalization including flavor mixings, and derive a evolution equation of the neutrino 
density matrix. In the section 11111 we will show how the predicted abundances of the light 
elements are modified when the reheating temperature is in the MeV range, and discuss 
their implications. Finally we present our conclusion in the section IIVI 



II. NEUTRINO THERMALIZATION 



In this section, we illustrate the formulation needed to follow the neutrino therma 
tion process. The case without the neutrino oscillations is studied by Refs. 



21 



22, 



iza- 
2J. 



Although subjects of study are different from this paper, issues of the neutrino spectrum 
evolution using momentum dependent Boltzmann equations in the early universe are treated 
in Refs. jl^l Q, 13, 0, [3, 33, Q • Our formulation almost goes in parallel with the no- 
mixing case and we use some of the useful techniques discussed in those papers. However, 
there is a very important exception that neutrino distribution functions have to be extended 



to neutrino density matrices |30| in order to include oscillations. 

First of all, let us explain our assumptions on the reheating process and the neutrino 
oscillations. We refer to the massive particles which reheat the universe, or inflaton, as <fi 1 . 
We assume <fi only decays into photons (they in turn produce electrons and positrons and 
they all thermalize very quickly by the electromagnetic interaction). In other words, the 
branching ratios to neutrinos or hadrons are assumed to be negligible and neutrinos are 
produced exclusively via the electron-positron annihilation. Then <fi is characterized simply 
by its decay rate V. We parametrize it by the reheating temperature Tr which is defined as 

T = 3H(T R ), (1) 

where H is the expansion rate of the universe. Here, we use the Friedmann equation 
H 2 = ptot/3Mpi where the reduced Planck mass M Pi = 2.435 x 10 18 GeV. The total en- 
ergy density p tot which consists of the radiation including photons, electrons, and three 



1 Hereafter we call <j> as inflaton even if it is not responsible for inflation. 
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species of neutrinos, is expressed as p to t = (5 , * 7I " 2 /30)T^ where the relativistic degree of 
freedom g* = 43/4. This leads to 

Tl _„ / T R V _ x 



T = 3.26—^- = 2.03 —4- sec" 1 . (2) 
Mpj V MeV / 

It should be noted that, even if the neutrinos are not fully thermalized, we stick to Eq. (J2J) 
as the definition of Tr to avoid unnecessary confusion. 

We consider three active f 
neglected as in Refs. 



avors of neutrinos: u e , v a and v T . When the oscillations are 

nnn 

|21|,|22J,|23J, there are only two sets of variables required to describe 



the neutrino evolution. They are the distribution functions for v e and which have to be 
distinguished since they interact differently with electrons; v e interacts via both neutral and 
charged currents while and v T has only the former interaction. Since and v T interact 
with electrons identically, we do not need to solve for the distribution function of v T . It is 
same as On the contrary, when we include the oscillations among them, and v T 

also have to be distinguished because their oscillations between v e are known to be different. 
Namely, we need to consider general three-flavor oscillations which require 9 real variables 
to fully describe our issue. However, if # 13 is zero, a simplification to two-flavor oscillations 
is possible by using non-mixing mass eigenstates u' and v' T instead of and v T 2 . Then 
v' and v e are described by two-flavor oscillations and v' T decouples from the oscillations. v' T 
just interacts with via neutral current and should behave as v T (or z/ M ) in the no-mixing 
case. 

Under those assumptions, the variables necessary for simulating thermalization of oscil- 
lating neutrinos are: the inflaton energy density p^, the photon temperature T, the v e - v 'u, 
two-flavor neutrino density matrix p p and the v' r distribution function f u ' T (p)- P P and f u > T are 
functions of neutrino momentum p. p v i s defined by expectation value of the product of the 
creation and annihilation operators |30j ]: 

*}(p)o,(q)) = (2vr) 3 5^ (p - q) [p p } tj , {ij} = {e,p}, (3) 

where aj(p) is the annihilation operator for negative-helicity neutrino of flavor i with mo- 
mentum p. Readers should bear in mind that the density matrix p p is just an extension 



2 Similar simplification is shown to be useful to analyze the evolution of neutrino-antineutrino asymmetries 
by Refs. 000. 



5 



of the occupation number to the mixed neutrinos, and should not confuse with the energy 
density, to which we refer as p u , p^, etc. Each diagonal component of p p is the neutrino 
distribution of the corresponding flavor, while the off-diagonal ones represent more subtle 
information on the correlation. For anti-neutrinos we can similarly define the density matrix 

Pp- 

ftJCp^Cq)) = (27r) 3 5( 3 )(p-q)[p p l., {i,j} = {e,fi}, (4) 



where &*(p) is the annihilation operator for positive-helicity neutrino of flavor i with momen- 
tum p. However, unless the lepton asymmetry is very large, we do not have to distinguish 
neutrinos from anti-neutrinos. In this case they are related to each other as p p = p^. We 
next derive the differential equations which govern their evolutions. 

We use scale factor a as a time variable and later we use y = pa instead of a momentum 
Then the time evolution equation for the neutrino density matrix p p is 3(j] 

Ha^ = -i[Sl(p),p p ] + I oM (p). (5) 

The matrix Q(p) represents both the vacuum oscillations and the refractive term. Neglecting 
lepton asymmetry, which is usually as small as the baryon asymmetry, it is written as 

n(p) - nv(p) - S 4^e, (6) 

where the Fermi coupling constant Gp = 1.16637 x 10 -11 MeV -2 , W boson mass mw = 80 
GeV. In the ultrarelativistic limit, £lv(p) is given by 

<V(P) = ^ UM2uT i ( 7 ) 

where M 2 , the neutrino mass matrix, and U, the matrix which relates mass eigenstates and 
flavor eigenstates, are 

\ m\ J \ — sin 6*12 cos 9 12 J 

We use the solar neutrino oscillation experiment values for neutrino parameters: m\ — m\ = 
7.3 x 10~ 5 eV 2 and sin 2 9 12 = 0.315 [38] 3 . The second term in Eq. (JBJ) comes from the 



The most recent result |2f| gives a slightly higher value of the mass squared difference: m| — m\ = 
7.9tgj x 10- 5 eV 2 . However, we have confirmed that our results do not change for m| — m\ in the error 
range. 



non-local effect of the W^-exchange interactions, and E is the energy density matrix of the 
charged leptons: 



E 




f (7/60)tt 2 T 4 

loo 



(9) 



where p e r S ) is the energy density of electrons (positrons) and we have assumed that neither 
muons nor taus exist in the plasma. 

For the collision term J co n, we consider the processes v + e ± <-> v+e^ and v + i> «-> e~ + e + . 
In calculating the collision term, we take electrons to be massless and neglect processes of 



scattering among neutrinos as Refs. [21 , 

Iveueipi) = TfijT J gfr ^ 2E4 + P2 ~ Ps ~ 



2 5 G| [4(pi -p 2 )(P3 • P4) F LL {V .< ' . !/■".( 



/I) P (2) ,.(3) p (4)' 



+4(pi -p 4 )(p 2 ■P3)F RR [ 



,M p( 2 ) , y (3) P (4) 



)] 



Iveue (pi 1 



2E l J 2E 2 2E 3 2E A K J yF1 F2 Fi FA! 



x 2 5 G 2 [4(p! • p 4 )(p 2 • p 3 ) e^, i/ 3 >, e< 4 >) 



+4(pi -p 2 )(p 3 ■Pi)F RR ( 



,(1) p(2) ,.(3) ? (4) 



)] 



I^eg(pi) = 2~E\ j ^2^^ (27r) 4 5 {4) (Pl + P2 - P 3 ~ P4) 

x 2 5 G 2 [4(p! • p 4 )(p 2 • p 3 ) i^z/ 1 ), z^ 2 ), e®, eW) 

+4(pi • p 3 )(p 2 • p 4 ) e^, gW)] , 

where we define dp = d 3 p/(27r) 3 , £!j = p°, and 

F ab (z/«,e( 2 ),z/( 3 ),e( 4 )) = l[(l-p w )G a p P3 G' 6 (l-/ e (p 2 ))/ e (p4)+h.c. 



F ab (u^,e^,^,e^) 



- P Pl G a {l - p P3 )G b f e (p 2 ){l - feipi)) + h.c] 
[(1 - p Pl )G a p P3 G b (1 - /g(p 2 ))/g(p4) + h.c. 
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-p pl G a (l - p P3 )G b f s (p 2 )(l - hipd) + h.c] 
5 P2 )G , 6 / e (p 3 )/g(p 4 ) +h.c. 

P Pl G a P P2 G b (l - / e (p 3 ))(l - /g(P4)) + h - c -] 



(10) 



(11) 



(12) 



(13) 
(14) 
(15) 



with G L = di&g(g L , g L ) and G R = di&g(g R , g R ). Here, g L = g L - 1 = sin 2 9 W - § and 
g R = sin 2 6*^ where sin 2 6 W = 0.23120 is the weak-mixing angle. / e (/g) is the distribution 



function of electrons (positrons). Hereafter we take p p = pi and fg = f e - Note that these 
collision terms coincide with those found in Ref. [19( if oscillations are absent (i.e., if off- 
diagonal components in the density matrices are zero). 

We further approximate that electrons obey the Boltzmann distribution and their Pauli 
blocking factors are neglected. Namely, in F's, we replace as f e {p) — > exp(— p/T) and 
1 — f e (p) — > 1- Then the collision terms above are reduced to one- dimensional momentum 



22j and the reduced expressions become equal to the 



integration by the technique in Ref. 
ones in the reference 4 in the limit of the zero mixing angle. 

In practice, since p p is 2x2 Hermitian matrix, it is convenient to expand it using Pauli 
matrices. Namely, 

3 

v P„ (V) 



Pp = E P ^)?> ( 16 ) 



i=0 



where 





<70 = ' ' , <Tl = , ^2 = , ^3 = • (17) 





On the right hand side of Eq. p p ] and J co u, are expanded similarly. We solve 

for the evolution of Pq ~ P3 and the distributions of u e and v' are in turn derived by 
f Ue = (P + Ps)/2 and /„/ = (P — Pz)/2. The evolution equations are formally written as 

Ra dPM_ = _ i n l (y)+/ ( (y), (18) 
da 

where % runs from to 3, Qi = Tr([0, p p ]f7i) and U = Tr(/ co n ai), and we have changed the 
variable p to y. 

We need to solve for the evolution of z/ , too. To this end, it is most simple to obtain the 
time evolution of f u > T from the z/^-component of Eq. (jSJ) with no mixing (which is given by 
omitting the first term on the right hand side) because the interactions of z/ with e ± are 
same as those of v'„. 



The right hand side of Eq. (A16) in Ref. [2^ has to be multiplied by 2 to be a correct equation. Due to 
this error, the right hand side of Eq. (3) in Ref. j^] has to be multiplied by 2. The right hand side of 
Eq. (12) in Ref. [22| has to be multiplied by 8 since it had already contained a typo of factor 4. In this 
occasion, we correct a typo in the right hand side of Eq. (8) in Ref. [22^: it has to multiplied by 2 (so that 
it is same as Eq. (2) in Ref. |2l|. 



s 



For the evolution of p^ and T, the equations are almost same as those found in Ref. |22|. 
We just need modifications due to our use of scale factor a as a time variable and discrimi- 
nation of from v T . For p^, it is given by 
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da aH a 
The equation of the total energy-momentum conservation is 

dptot 3 



(19) 



da a 

where the total energy density and the total pressure are given by 

Ptot = P4> + P~i+ Pe± + Pu e + Pul + Pu' T 

vr 2 T 4 2 r , 2 E P 



(ptot + Pot), (20) 



+ —rr- + s / d PP 



15 vr 2 J exp(P e /T) + 1 



1 f°° 

-5-4 / dyy 3 (f„ e + f„,+U), (21) 

7T a Jo 



3 tot EE P 7 + P e± + P„ e + P„, + P„, , 

vr 2 T 4 2 r 00 , p 4 
+ ^5 dp 



45 3vr 2 i (exp(P e /T) + 1) 

^ / dyifiU + fy + U) (22) 



3tt 



with the electron energy E e ee \Jm1 + p 2 . The evolution equation for T is derived from 
Eq. flU}: 



da V^T+^T-J \a^ + a (Pe±+Pe±; off"" 

"'"^a 4 j o ^ ^ ^\ da ^ da ^ da^j\ ^ ^ 

Finally, the expansion rate is 

H ee 1^ = (24) 

To integrate the differential equations, since the equations for f v (y) are stiff, we used semi- 
implicit extrapolation method [40]. Using the Ref. |40j's implementation which incorporates 
adaptive stepsize control routine, we were able to evolve the neutrino density matrices very 
efficiently. We followed the evolution well after the electron-positron annihilation ends and 
fuivYs become constant. 



As for the initial condition, we have to make the inflaton energy density dominate the 
universe at first. As long as p^ is much larger than radiation energy density (~ T 4 ), evolution 
afterward does not depend on their precise values. In this paper, we adopt rather realistic 
relation between p^ and p ra( j, 

Prad = ^rM PlP y 2 , (25) 

which derived from the analytic solutions during the epoch of coherent oscillations 



III. RESULTS AND COSMOLOGICAL IMPLICATIONS 

In this section, we present the results of our numerical calculation for neutrino thermal- 
ization and consider its implications for cosmology. We evolve the neutrino density matrices 
with various values of the reheating temperature Tr and investigate how the neutrino dis- 
tribution functions, neutrino energy densities and big bang nucleosynthesis depend on T R . 



Along with the neutrino thermalization with osci 



lations which have been studied in Refs 



21 



22, 



lations, we show the results without oscil- 
2311 and elucidate the neutrino oscillation 



effects on a low reheating temperature scenario. Our results when the oscillations are omit- 
ted turn out to be consistent with those of previous papers. We find that the inclusion of the 
oscillations most characteristically alters the 4 He synthesis and its abundance varies with 
respect to T# quite differently from the no oscillation case . 

A. Neutrino distribution functions 

We show the final neutrino distribution functions in Figs. ^ (a)-(d) for the cases of Tr = 
15 MeV and 2.5 MeV respectively with and without the oscillations. We see from Figs. ^ 
(a) and (b) that when the reheating temperature is sufficiently high, all the neutrino species 
are thermalized regardless of the oscillations. 

However, for the case of lower reheating temperature, the oscillations significantly matter 
as seen from comparing Figs. (c) and (d) : f Ue and f u r are almost equalized by the solar 
mixing. When the oscillations are neglected, f Ue becomes much larger than f v as shown 
in Fig. [T] (c) because v e is produced by the charged current interaction in addition to the 



neutral current interaction but and u T are produced only by the latter 
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FIG. 1: The final distribution functions of neutrinos, (a) and (c) are cases for no oscillations 
[v e is displayed by solid lines and by dashed lines) and (b) and (d) incorporate the oscillations 
[y e is displayed by solid lines, u' by dashed lines and v' T by dot-dashed lines). The equilibrium 
distributions are drawn by dotted lines in order to show how much they are thermalized. For 
Tr = 15 MeV, in (a) and (b), whether the oscillations are present or not, all the lines overlap and 
this means every neutrino species is fully thermalized for high reheating temperature. For Tr = 2.5 
MeV, in (c) and (d), distributions are away from equilibrium form. When the oscillations are taken 
into account, distributions of v e and v' get close as seen in (d). 



there are the flavor mixings, v e and v' can convert into each other, u' is now produced 
also by the oscillations from u e which exists more than u' so f v i increases compared to no 
oscillation case. On the contrary, f Ve becomes smaller when the oscillations are included 
naturally because v e oscillates into v'. However, this deficit is to some extent filled by the 
v e production from the thermal plasma so the neutrinos are produced more in total under 
the existence of the oscillations. This is seen in Fig. |21 which shows clearly that f Ve + f v > 
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FIG. 2: We draw the sums of the distribution functions, f„ + /„ (no oscillation) and /„ + f v i 
(including oscillation) with the dashed line and the solid line respectively. The latter is larger 
showing that the oscillations make the thermalization more efficient in total. 

(including oscillation) > f Ve + f v (no oscillation). 
B. Effective number of neutrinos 

Let us discuss our results of the neutrino thermalization in terms of neutrino energy 
density. This is often expressed using the effective number of neutrinos N v . This number is 
observationally relevant to the CMB power spectrum and large scale structure. It is given 



where the summation is taken for v — u e , and v T when the oscillations are not included 
and v = u e , u' , and u' T when we consider the oscillations. We define p^std using the photon 
temperature T as 



by 



N v = 



(26) 




(27) 
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FIG. 3: The effective neutrino number N v as a function of the reheating temperature Tr (shown 
on the bottom abscissa) or the decay width V (shown on the top abscissa). The cases with and 
without the oscillations are drawn respectively by the solid and dashed lines. The horizontal line 
denotes N v = 3.04 with which N v for high Tr should coincide (see the text). 



which corresponds to the neutrino energy density assuming that neutrinos are completely 
decoupled from the rest of the thermal plasma before the electron-positron annihilation takes 
place. If this assumption is exact, N u would be 3. It is actually a very good assumption but 
detailed calculations on the entropy transfer from electrons to neutrinos have shown that 
p„'s are slightly larger than p^std and N v = 3.04 QQHQ- 

We calculate p v by integrating the final neutrino distribution functions such as presented 
in Fig. ^ and derive N v as a function of the reheating temperature Tr. The result is shown 
in Fig. El For Tr > 10 MeV, N v asymptotes the value 3.04 which indicates thermalized 
neutrino distributions. This is regardless of the neutrino oscillations and consistent with 
Fig. [T] (a) and (b) discussed in section 1111 Al For the smaller values of T R , the inclusion of 
the oscillations make N u larger as expected from Fig. |21 This effect is most conspicuous for 
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T R = 2 ~ 5 MeV and changes N v up to ~ 0.2. 

Fig. El enables us to constrain Tr by using the limits on the effective number of neutrino 
species from cosmological observations such as CMB and galaxy surveys. Recent papars, 



Refs. 



41 



42 



43 



44j . derive the lower limit to be 0.9 ~ 1.9 (these are the limits obtained 



without resorting to BBN. Some of them have also reported more stringent limits obtained 
using data combined with observed Y p . However, since they assume Fermi-Dirac distribution 
for neutrinos and only modify the Friedmann equation when they calculate Y p , we cannot use 
those limits. This point is discussed in section ITlI CI in detail). If N u > 0.9 is adopted, the 
bound on the reheating temperature is Tr > 1.69 MeV with the oscillations and Tr > 1.74 
MeV for no oscillation case. 



C. Light element abundances 

We now investigate how the big bang nucleosynthesis is affected by the non-thermal neu- 
trino distributions and/or the neutrino oscillations. We calculate the light element (D, 4 He 
and 7 Li ) abundances as functions of Tr, again with and without the neutrino oscillations. 
The cosmological effects of incomplete neutrino thermalization is most strikingly seen in 
4 He abundance since electron-type neutrinos play a special role in determining the rate of 
neutron-proton conversion during BBN. This has been already known from the previous 



papers Refs. |21|,|22| in which the oscillations are neglected, but we find that the neutrino 
oscillations prominently matter in regard to the T^-dependence of 4 He abundance. 

We show how Y p varies with respect to Tr in Fig. 0] This is calculated by plugging 
ihe solutions of the evolution equations derived in section III] into the Kawano BBN code 

"i , n 

45| (with updated reaction rates compiled by Angulo et al. |46j). Required modifications 
are the temperature dependence of the neutron-proton conversion rates, T n ^ p and T p ^ n , 
and the evolution equation for the photon temperature. The calculation of r n ^ p (see e.g. 
Ref. (4?]]) involves the integration of the electron neutrino distribution function f Ve which 
does not necessarily take the Fermi distribution form in our case. For the photon temperature 
evolution, the contributions from and neutrinos are supplemented in the same way as 
Eq. 03). 

There are two effects caused by incomplete thermalization of neutrinos competing to 
make up the dependence of Y p on Tr as shown in Fig. EJ slowing down of the expansion rate 
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FIG. 4: The He abundance (mass fraction) Y p as a function of the reheating temperature Tr 
(shown on the bottom abscissa) or the decay width T (shown on the top abscissa). The cases with 
and without the oscillations are drawn respectively by the solid and dashed curves. Thinner curves 
are calculated with fermi distributed neutrinos with N v of Fig. |3] (namely, only the change in the 
expansion rate due to the incomplete thermalization is taken into account). The horizontal line 
represents "standard " Y p calculated by BBN with neutrinos obeying the fermi distribution and 
N u = 3.04. The baryon-to-photon ratio is fixed at rj = 5 x 1CP 10 . 

and decreasing in r n ^ p . The former is just a result of the decrease in the neutrino energy 
density (of all species). The latter is due to the deficit in f Ue . They compete in a sense that 
they work in opposite way to determine the epoch of neutron-to-proton ratio freeze-out: the 
former makes it later and the latter makes it earlier. Then, the competition fixes the n-p 
ratio at the beginning of nucleosynthesis and eventually determine Y p . Roughly speaking, 
for larger Tr, the former dominates to decreases Y p but, for smaller Tr, the latter dominates 
and increase Y p . This is clearly seen in the case without the oscillations but not for the 
case including the oscillations because the incompleteness in the v e thermalization is made 
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T (MeV) 

FIG. 5: The weak interaction rate T n ^ p and the expansion rate H as functions of temperature, 
where r ra __> p and H first become equal. We plot for Tr = 2.5 MeV with and without the oscillations. 
For Tr = 15 MeV, the oscillations do not make any difference. 

severer by the mixing (see the panels (c) and (d) in Fig. ^) and this effect dominates already 
at high T R . 

Before going forward, it may be worthwhile to look slightly more into the explanation of 
the T R - dependence of Y p . First, let us forget about modifying T n ^ p or temperature evolution 
and just calculate 4 He abundance using thermally distributed neutrinos with N u : s indicated 
in Fig. El for each value of Tr. This corresponds to including the effect of slowing down the 
expansion rate due to the incomplete thermalization but neglecting the electron neutrino 
deficiency. Accordingly, lowering Tr only acts to delay the n-p ratio freeze-out and decrease 
Y p (shown by the thinner curves in Fig. 0J). In actual low- reheating temperature scenario, 
a lack of v e reduces T n ^ p . This counterbalances the effect of slowing-down expansion and 
boosts Y p in total at lower Tr. To see this is really the case, we plot T n _> p for some values of 
Tr in Fig. We see that T n ^ p is smaller for lower Tr which is attributed to less thermalized 
u e . It is also instructive to calculate the neutron-to-proton ratio freeze-out temperature T np , 
which we define by r n ^ p (T np ) = H(T np ), to confirm where the competition settles. This 
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FIG. 6: (a) T np , freeze-out temperature of the neutron-to-proton ratio, and (b) 2/[l + 
exp(Am/T np )], as functions of Tr. 



is shown in Fig. El (a) and we see that at low Tr, the decrease in T n ^ p wins to make 
T np higher (in the case with the oscillations, this seems to win for every Tr and T np rises 
monotonically as T R decreases). We note that the figure well reproduces the profile found 
in Fig. 0] This resemblance becomes more meaningful by plotting instead the quantity 
2/[l + [n p /n n )f] = 2/[l + exp(Am/T np )], the usual estimation of 4 He abundance from the 
neutron-to-proton ratio at the freeze-out value, which is shown in Fig. El (b). Although the 
figure is not exactly same as Fig.lUbecause free decays of neutrons are not considered, we see 
that the 1^,'s dependence on Tr is sufficiently understood from this estimation. When the 
neutron free decay is properly taken into account, the estimation for Y p decreases from the 
values indicated in Fig. El (b) . For lower Tr, since the time between T np and the start of the 
nucleosynthesis (T w 0.07 MeV) is longer (this in turn is explained by the smaller expansion 
rate due to less neutrino energy densities), this decrease should be larger. Therefore, on 
including the neutron free decay, Fig. El (b) would be tilted toward left (smaller Tr) side 
and should look more like Fig. |3J In particular, the minimum found for the case without 
the oscillations should be located at lower T R when the free decay is included. 

We have so far discussed the 4 He synthesis features common to the 1ow-Tr universe with 
and without the neutrino oscillations, but we would rather like to emphasize that there is 
a striking difference between them. This is most clearly visible in Fig. EJ when we include 
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the oscillations, Y p does not decrease if we lower Tr. This is somewhat surprising because, 
at the same time, N v becomes smaller (see Fig. |3J). This means that, in the case with the 
oscillations, the effect of slowing down cosmic expansion (as represented by decreasing N u ) 
is completely overcome by the decrease in T n ^ p for all Tr. The reason why this happens is 
that since the oscillations convert electron neutrinos into muon neutrinos, the deficiency in 
electron neutrinos is made severer (see Fig- 13) - Moreover, why this matters for 4 He synthesis 
is that it is exclusively sensitive to the u e distribution function which determines T n ^ p . On 
the other hand, the structure formation is affected only by the energy density so it does not 
distinguish neutrino flavors. Since only their sum matters, the oscillations scarcely make 
difference (see Fig. |HJ). Therefore, BBN, especially when the neutrino oscillations are taken 
into account, turns out to be unique probe of low reheating temperature scenario. Next, we 
proceed to compare the predictions of the scenario with the observed abundances. 

On comparing the predictions of low reheating temperature scenario with the observed 
abundances, we need to vary the baryon-to-photon ratio, i], which is the input parameter 
for the standard BBN calculation, in addition to Tr. In Fig. we show contour plots for 
abundances of light elements, D, 4 He and 7 Li, against 77 and Tr. Since contours tend to 
be parallel to each other, we see that how abundances vary with respect to Tr has little 
dependence on 77. In particular, for 4 He, features found in Fig. H] seem to appear at every 77. 
We notice, in Fig. [7| (c) and (d), that the oscillations almost do not make difference for D 
and 7 Li abundances. In the figure, we also indicate observed values taken from Ref. for 
4 He, from Ref. [3] for D, and from Ref. for 7 Li: 

Y p = 0.238 ±0.002 ±0.005, (28) 
D/H = 2.78± 1 ; 4 ^ x 10~ 5 , (29) 
'Li/H = 1.23^;1 x !0~ 10 (95%) (30) 



7- 



In Eq. (|28j). the first uncertainty is statistical and the second one is systematic. Their root- 
mean-square, [(stat.) 2 + (syst.) 2 ] 1 / 2 , is adopted as overall la error. In this paper, we do not 
consider the 7 Li data since its systematic error is under debate at present, but show it just 
for reference 5 . 



5 It is known that the baryon density derived from Eq. I|28(l is somewhat lower than one from Eq. i|29|l . It 
was widely believed that N v < 3 decreases Y p and ameliorates this tension (see e.g., Ref. However 
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We immediately realize from Fig. [7| (a), (b), (c) that inclusion of the oscillations leaves 
smaller room for the low reheating temperature scenario. In other words, the parameter 
region allowed from D and 4 He measurements is smaller for the case with the neutrino 
oscillation. We can see it more clearly by x 2 -analysis whose results are shown in Fig. |SJ The 
lower bound on Tr at 95% confidence level in the tj-Tr plane is 1 MeV for the case of no 
oscillations but tightened to be 2 MeV for the case incorporating the oscillations 6 . 

IV. CONCLUSION 

In this paper we have investigated the MeV-scale reheating scenario wherein the thermal- 
ization of neutrinos could be insufficient. We have paid particular attention to the oscillation 
effects on the thermalization processes of neutrinos, and solved numerically the momentum 
dependent Boltzmann equations for neutrino density matrix, fully taking account of neutrino 
oscillations. In contrast to the widespread picture, we have found that 4 He abundance does 
increase while the effective neutrino number N v decreases. The reason is simple; the neutrino 
oscillations reduce the number density of v e , due to which the neutron-proton transforma- 
tion decouples earlier. This effect cancels and even overcomes that of the decrease in the 
expansion rate; only the latter effect has been usually taken into account when discussing 
the effect of N v on the light-element abundances. Therefore we would like to stress that it 
is indispensable to take into consideration the oscillation effects, to set a lower bound on 
the reheating temperature by using the BBN. As a reference value, we quote our results; 
Trh ^ 2 MeV or equivalently N v > 1.2 obtained by using the observational data on the 4 He 
and D abundances. 

What are then the distinct predictions of the MeV-scale reheating? Clearly, they are: 
both larger Y p and smaller N v compared to their standard values; if both the observed Y p 
and N v suggest the same Tr by the relations shown in Figs. El and IH they would serve as 

we now know this is not feasible simply by lowering the reheating temperature if we properly take into 
account neutrino oscillations. 
6 Recently, analysis of the 4 He abundance by Ref. [52( suggests Y p = 0.249 ±0.009 |53|. This is higher than 
the value of Eq. I|28|) mainly due to the different treatments of stellar absorption. Although, at present, 
such large uncertainty does not allow us to derive any meaningful lower bound on Tr. However, higher 
Y p is interesting for MeV-scale reheating scenario. Should future research yield Y p > 0.25, Tr ~ O(MeV) 
would be favored. 
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(a) Yp (Inc. oscillation) (b) Yp (No oscillation) 




FIG. 7: Contour plots for the light element abundances. 4 He mass fraction is plotted in (a) with 
the oscillations and in (b) without. D and 7 Li are plotted respectively in (c) and (d) where dotted 
lines express the case with the oscillations and dashed lines express the case without. Shaded areas 
represent uncertainties in the observed abundances expressed in eqs. (j28j) ~ (|3Uj) (for D and 7 Li, 
they are drawn against the contours considering the oscillations). Darker areas are for 1 a and 
lighter for 2 a. 

decisive evidence for the MeV-scale reheating 1 . 

At last, let us comment on the validity and possible extension of the present work. As ex- 
plained in section m we have neglected the self-interactions of neutrinos. Such simplification 
is considered to be valid due to the following reason. Since self-interactions cannot change 

7 According to Ref. , we can determine both Y p and N v with future CMB observations such as Planck. 
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FIG. 8: x 2 contour plot using data of D and 4 He. For no oscillation case, the allowed regions at 
and 95% confidence levels are drawn with solid and dashed lines. For the case with oscillations, 
the 68% allowed region does not appear and only the 95% region is indicated by the shaded area. 



the total energy stored in the neutrino sector, they affect only the momentum distribution 
of neutrinos. On the other hand, it should be noted that we have taken into consideration 
the neutrino-electron (ve) scattering, which also shifts the neutrino momentum distribution 
toward kinetic equilibrium at the rate of the same order of magnitude as the vv scatter- 
ing. However we have checked that our results do not change at all even if we increase 
the ve scattering rate a few times larger than the standard one. Considering that the vv 
scattering rate is further suppressed due to the deficit in the neutrino number, we are sure 
that the self-interactions have only a minor effect in the neutrino momentum distribution. 
Still, the self- interactions have a potential effect on the number density of v e through e.g., 
v e v e «-> v^r) Vftfr) ■ Furthermore, nonzero #13 can have a similar effect; in this case it is neces- 
sary to perform three generation analysis. Nevertheless we believe that our main conclusion 
is robust, since these extensions, too, are expected to decrease the number density of v e , 
further increasing the 4 He abundance. Of course the quantitative improvement should be 
necessary and the full analysis on these points will be presented elsewhere j^. 
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